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I. INTRODUCTION 

One of the most fundamental properties of a fluid is the surface tension at the liquid- vapor 
interface. It would seem that such a fundamental property would be an ideal candidate for 
study via computer simulation. However, the determination of the surface tension from 
simulation turns out to be frought with difficulties so that even today there is still a sub- 
stantial amount of effort directed towards the development of more reliable algorithms and 
the refinement of the reported values even for the paradigmatic case of a simple fiuid mod- 
eled with the Lennard- Jones interaction pj. One of the primary difficulties is that in all 
simulations, the potential is truncated at some finite range and it happens that the surface 
tension is very sensitive to the value of the cutoff. For that reason, an important part of 
the development of algorithms has focused on the calculation of the corrections needed to 

Bit the infinite-ranged limit from data obtained using a truncated potential (see, e.g., ref. 
, y]). This sensitivity is therefore a nuisance when the goal is to get the infinite range 
result, but in other ways it can made useful. In particular, one of the important reasons to 
determine the surface tension from simulation is that it provides a baseline against which 
theories of inhomogeneous liquids can be tested [3|,|j]. For this application, the sensitivity of 
the surface tension to the range of the potential can be used as a test of the generality of a 
theory which was probably motivated in the first place by its agreement with some existing 
simulation data. Furthermore, there has recently been a significant increase in interest in 
short-ranged potentials in their own right. This is due to the fact that certain complex fiuids, 
in particular globular proteins, can, in a first approximation, be modeled as a simple fiuid 
with a very short ranged interaction |5|. It is therefore interesting to study the properties 
of fluids with these kinds of interactions and to test that existing theories work in this new 
domain of interest. For these reasons, we present in this paper a systematic study of the 
dependence of the surface tension of a Lennard- Jones fiuid as a function of the range of the 
potential. 

In this paper, we describe the results of Monte Carlo (MC) simulations of a Lennard- Jones 
fiuid with the potential truncated at several different points. We have chosen to truncate 
and shift the Lennard- Jones potential, VLj{r), so that the potential used in this work is 
f (r; Tc) = VLj{r) — Vi,j{rc) for r < Tc and f (r) = for r > r^. If this shift is not performed, 
then there is an impulsive contribution to the pressure when atoms move across the r = Tc 



boundary that would have to be taken into account 6|. We do not shift the force, i.e. we 
do not use w(r; Tc) = VLj{r) — Vlj{j-c) — (r — rc)v'^j{rc) with f^j(r) = dvLj{r)/dr inside the 
cutoff, as is usually done in molecular dynamics simulations to avoid impulsive forces: our 
potential is truncated and shifted but the force is not shifted. This choice was made in order 
to allow for comparison with previous MC studies. 

In the simulations a slab of liquid is bounded on both sides by vapor. The surface tension 
is determined using the method of Bennett [g, [3] as there seems to be some evidence that 
this method is more robust than other commonly used techniques |8(]. It is often the case 
that the quantity of interest is the surface tension for the infinite-ranged potential. Since 
simulations almost always make use of truncated potentials, various techniques have been 
developed to approximate the so-called long range corrections, i.e. the difference between 
quantities calculated with the truncated potential and the infinite-ranged quantities [2]. We 
do not include any such corrections here since our goal is actually to study the truncated 
potentials. Thus, each value of the cutoff defines a different potential with its own coexistence 
curve and thermodynamics. 

In Section II, we present the simulation techniques used in our work. Section III con- 
tains a discussion of our results including a comparison to previous work. Since one of the 
motivations for this work is to provide a baseline for testing theories of the liquid state, we 
illustrate this by comparing our results to Density Functional Theory (DFT) calculations 
and by testing the law of corresponding states. We give our conclusions in the last Section. 



II. SIMULATION METHODS 

Simulations are performed with a standard Metropolis Monte-Carlo algorithm (MC-NVT) 
for a system of A^ = 2000 particles of mass m at temperature T in a volume V = L^LyL^ 
where L^,, Ly, and Lz are the dimensions of the rectangular simulation cell. Periodic bound- 
ary conditions are used in all directions. Particles interact via the Lennard- Jones potential, 

-W-<(7)'^-(7)°) (^) 

which is truncated and shifted so that the potential simulated is 

v[r) = <( (2) 

: r > Tr 



where Tc is the cutoff radius. Each simulation starts from a rectangular box (L^. = Ly = 
L, Lz = 4:L) filled with a dense disordered arrangement of particles (p* = pa^ = 0.8) 
surrounded along the z-direction by two similar rectangular boxes containing particles in a 
low density state (p* ~ 0.01) . The total simulation box has sides of length L^ = Ly = 9.15(7 
and Lz = 109.63(7. The liquid film located in the middle of the box has a thickness Az ~ 27a 
so that the two interfaces do not influence each others. The system is first equilibrated during 
5 X 10^ Monte-Carlo cycles (one cycle = N updates) after which the positions of the particles 
are saved every 20 cycles during 5 x 10^ cycles. This ensemble of 2.5 x 10^ configurations is 
used to compute the density profile and the surface tension by the Bennett's method. 

Although several methods are available for the computation of the surface tension, the 
Bennett's approach has been chosen because of its accuracy[8]. In the Bennett's method the 
calculation of the surface tension follows from the definition 

\(^^J N,V,T 
where F is the free energy and A is the area of the liquid-vapor interface. In its imple- 
mentation the method requires that one performs two simulations: one for system of 
interface area Aq, and another for system 1 of interface area Ai = Aq + A A. In this work 
AA/A = 5 X 10~^. The free energy difference AF between the two systems is evaluated by 
the method of acceptation ratio which starts with the computation of Aii^oi = -^oi ~ -^oo 
which is the difference between Eqq, the energy of a configuration of system 0, and Eqi, the 
energy of a new configuration obtained from the previous one by rescaling the positions of 
the particles (a, l9|] : x' = x {Ai/Ao)2^ y' = y (Ai/Ao)2^ and z' = z{Ao/Ai). Similarly one 
computes AEiq = Eiq — En obtained from a configuration of system 1 following an inverse 
rescaling of the positions. AF is obtained by requiring that 

Y, fi^Eoi ~AF) =J2 fi^E.o + AF) (4) 

no ni 

where ^^o (^m) ^^ ^ ^^^ °^^^ ^^^ configurations of systems (1), and f{x) = (1 + 
exp(/9x))^^. Then, taking into account the fact that the system contains two fiat interfaces, 
the value of the surface tension is given by 7 = AF/{2AA) . 




FIG. 1: (Color online)The surface tension as a function of temperature for two different cutoffs. 



The open circles are our data, the filled circles are from Duque et al 
Mecke et al[l[, the diamonds are from Potoff et al 



ll|, the squares are from 



12] and the triangles are from [10]). Note that 



the Mecke and Potoff data both include long-ranged corrections. 



III. RESULTS 



A. Comparison to previous results 



Our results for the surface tension as a function of the cutoff are given in Table 1. Note 
that all quantities are reported in reduced units so that the reduced temperature is T* = T/e, 
the reduced cutoffs are r* = rda and the reduced surface tension is 7* = 7cr^/e. In Fig. [1] 
we show our results for cutoffs of r* = 2.5 and 6.0 compared to the MC data of Have and 

n n 

Bruin [lCf| for the shorter cutoff and to the MD data of Duque et al[ll| (who appear to shift 



the forces) and Potoff et al[l2i] and Mecke et al[l|. The latter two are shown even though 
they include long-ranged corrections. Our data are seen to be very consistent with the MC 
data obtained without long-ranged corrections and to lie slightly below the corrected data, 
as expected. 



B. Comparison to DFT 

In Fig. [21 we compare our results to the predictions of a recently proposed Density 
Functional Theory modelj^. The DFT requires knowledge of the bulk equation of state and 
the figure shows results using two different inputs: the 33-parameter equation of state of 
Johnson, Zollweg and Gubbins (JZG)jl3| and first order Barker-Henderson thermodynamic 
perturbation theory [ij, [l5|. Both versions of the theory are in good qualitative agreement 
with the data, showing the decrease in surface tension as the range of the potential decreases. 
It might be thought that use of an empirical equation of state should automatically give 
superior results to an approximation, like thermodynamic perturbation theory, but this is 
not necessarily the case since the equation of state is fitted to data for the infinite-ranged 
potential. The finite cutoff is accounted for using simple mean-field correct ions J6|, [l3| and 
these become increasing inaccurate as the cutoff becomes shorter and, for fixed cutoff, as the 
fluid density becomes higher. The latter condition means, in the present context, increasing 
inaccuracy as the temperature decreases. Both of these trends are confirmed by the figure. 
The decrease in accuracy with decreasing cutoff can be seen in the fact that the critical 
point (corresponding to the temperature at which the surface tension extrapolates to zero) 
is less accurately estimated for the smaller cutoffs than for the larger cutoffs. 

The perturbation theory, on the other hand, takes the cutoff into account more accurately 
and consistently so that no strong change in accuracy is expected as the cutoff decreases. 
However, the theory itself is expected to be less accurate for higher densities so again a 
drop in accuracy with decreasing temperature would be expected and that is indeed seen in 
the figure. Furthermore, perturbation theory is in general going to be inaccurate near the 
critical point as it does not take into account renormalization effects which tend to lower 
the critical point. These effects are less pronounced for shorter-ranged potentials and indeed 
the perturbation theory seems more consistent with shorter-ranged potential. 

C. Corresponding states 

The principle of corresponding states is a generalization of the results of the van der 



Waals equation of state 16j. The idea is that the properties of simple liquids should be 



universal functions of the state variables, density and temperature, scaled to the critical 




FIG. 2: (Color online) Comparison of our simulation data to a DFT model[4]- The panel on the left 
shows the results of the theory using an empirical equation of state while the results on the right 
were obtained using thermodynamic perturbation theory.The lines are ordered from the smallest 
cutoff (lowest lines) to the largest cutoff (highest lines) and were calculated for r* = 2.5, 3, 4, 6, oo. 
The data is represented by circles(2.5), squares(3), diamonds(4), filled diamonds (4 - larger system) 
and triangles(6). 



point. In this section, we test this hypothesis by applying it to the surface tension. The 
first step is therefore to determine the critical temperatures and densities of the various 
truncated potentials. Since the theoretical calculations require as input an equation of 
state, the critical points are easily determined. To determine them from the simulations, 
we took five independent averages over 5000 configurations and fitted the density profiles 
in each case to a hyperbolic tangent and then from these extract the coexisting vapor and 
liquid densities at each temperature. The five values obtained at each temperature were 
averaged and the variance used as an estimate of the errors in the values. The critical 
temperature was then estimated b y u sing the lowest order renormalization group (RG) 



result, i(p; 



Pv) 



A{T, - T) 



0.325 



17 



18| . In some applications 19j , higher order terms 



are included but we did not feel that the accuracy of our data warranted use of anything 



TABLE I: Surface tension determined from simulation as a function of temperature for different 
cutoffs. 



Temperature 


r* = 2.5" 


r* = 3" 


r*=4" 


r*=4^ 


r* = 6^ 


0.70 


0.584 (27) 


0.770 (21) 


0.964(46) 


0.914(30) 


1.070(13) 


0.72 


0.561 (26) 


0.726 (25) 


0.899(22) 




1.034(7) 


0.75 


0.511 (20) 


0.698 (28) 


0.825(31) 




0.959(8) 


0.80 


0.421 (19) 


0.608 (16) 


0.748(25) 


0.736 (32) 


0.847(13) 


0.85 


0.315 (13) 


0.480 (12) 


0.633(24) 






0.90 


0.228 (13) 


0.384 (15) 


0.542(16) 


0.484(33) 


0.638(8) 


0.95 


0.181 (11) 


0.314 (12) 


0.443(22) 






1.0 


0.106 (7) 


0.234 (8) 


0.348(11) 


0.313 (59) 


0.438(15) 


1.05 




0.156 (11) 


0.258(23) 






1.1 




0.111 (16) 


0.200 (12) 




0.261(13) 


1.15 




0.074 (10) 


0.108 (14) 






1.2 




0.054 (6) 


0.067 (15) 


0.063 (19) 


0.105(5) 



"using approximately 2000 atoms, 
^using approximately 8000 atoms. 



but the lowest order function. The critical density was then estimated using the law of the 
rectilinear diameter, j{pi + Pv) = Pc + -S(Tc — T) [16]. There are again higher order corrections 
to this formula which can be calculated using RG methods, but for the reasons just given, 
we have not attempted to include them. The results of these fits are illustrated in Fig. [3] 
and summarized in Table 2. The largest errors in this procedure are in the determination 
of the critical density. 

Figure m shows the surface tensions, as determined from simulation and theory using the 
JZG equation of state, scaled to the critical density and temperature as a function of distance 
from the critical temperature. Despite the wide range of cutoffs and the mixture of data 
from simulations and theory, it is nevertheless seen that the data do in fact obey the law of 
corresponding states to a good approximation. However, the same scaling of the theoretical 
calculations using the equation of state from thermodynamic perturbation theory, shown in 
Figure El does not give a single curve. While the data for the shorter cutoffs, r* = 2.5 and 
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T 



FIG. 3: (Color online) The fit of the difference in liquid and vapor densities, as determined from 
simulation (symbols), to the RG functional form (lines). The data are shown as circles {Re = 6.0), 
open squares {Re = 4.0, larger cell), filled squares {Re = 4.0, smaller cell), diamonds {Re = 3.0) 
and triangles {Re = 2.5). 

r* = 3 appear to coincide, the data for the larger cutoffs does not. This appears to be due, 
at least in part, to the fact that the estimate of the critical density as a function of the 
cutoff calculated using the perturbation theory is not monotonic (see Table [IT]) which is at 
odds with the quantities as determined from simulation which clearly are monotonic in the 
cutoff. 



IV. CONCLUSIONS 



We have presented our determination of the liquid-vapor surface tension in a Lennard- 
Jones fluid as a function of the range of the potential. The data give a systematic picture of 
the variation of surface tension with the cutoff and are in agreement with previous studies. 
It is hoped that this can serve as a useful benchmark for the development of theories of 
inhomogeneous liquids. Indeed, the results were compared here to calculations made using a 



TABLE II: The critical points for the LJ potential truncated at different values. The theoretical 

n 

values were determined using the empirical JZG equation of state|13f| (JZG) and the first order 
Barker-Henderson perturbation theory (BH). 



Rl 


MC 


Theory - 


JZG 


Theory - 


BH 




T* 


P*c 


T* 


Pc 


t: 


Pc 


2.5 


1.10 (1)" 


0.31 (9)" 


1.04 


0.26 


1.18 


0.325 


3.0 


1.18 (l)"^ 


0.31 (7)'^ 


1.15 


0.28 


1.27 


0.342 


4.0 


1.26 (1)'^ 


0.31 (5)^^ 


1.25 


0.32 


1.34 


0.341 


4.0 


1.25 (2)^ 


0.30 (6)'' 


— 


— 


— 


— 


6.0 


1.30 (2)^ 


0.32 (9)'' 


1.29 


0.35 


1.38 


0.341 


oo 


1.31 " 


0.317 '^ 


1.311 


0.351 


1.40 


0.312 



"using approximately 2000 atoms, 
''using apprpxlmately 8000 atoms. 
"^Fromref. [20 1. 



recently developed DFT and the strengths and weaknesses of the theory are evident: while 
it gives a good semi-quantitative estimate of the surface tension for all cutoffs, errors on the 
order of 10% are present indicating that further improvement is possible. 

We have also tested the law of corresponding states by showing our results from both 
simulation and theory scaled to the critical density and temperature. For the simulation 
data and the theoretical calculations based on an empirical equation of state, the law of 
corresponding states appears to be obeyed. However, the calculations based on the equa- 
tion of state from first order perturbation theory do not appear to scale well at all. This 
failure appears to be due to poor behavior of the critical density as a function of the cutoff 
and indicates that care must be exercised before using the law of corresponding states to 
extrapolate calculations. 
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FIG. 4: (Color online) The scaled surface tension, 7** 



Tcp: 



2/3 



, as a function of distance from the 



critical temperature. The left panel includes the theoretical curves, shown as full line {R* = cxd), 
dotted hue (i?* = 8), dashed hue {R* = 6), dash-dot line (i?* = 4), dash-dot-dot hue {R* = 3), 
and line+circles (R* = 2.5), and the simulation data, shown as circles {R* = 6), filled squares 
{R* = 4, 2000 atoms), open squares {R* = 4, 8000 atoms), diamonds {R* = 3) and triangles 
{R* = 2.5). The right hand panel shows only the data from simulation as well as the estimated 
error. In both cases, the thin line is a best fit to all of the data (theory and simulation) of the form 
7* = 7o*(l - T/T,) with 7o* = 3.41. 

(ARC 2004-09). 



[1] M. Mecke, J. Winkelmann, and J. Fischer, J. Chem. Phys. 107, 9264 (1997). 

[2] D. Duque and L. F. Vega, J. Chem. Phys. 121, 8611 (2004). 

[3] K. Katsov and J. D. Weeks, J. Phys. Chem. B 106, 8429 (2002). 

[4] J. F. Lutsko, J. Chem. Phys. 128, 184711 (2008). 

[5] P. R. ten Wolde and D. Frenkel, Science 77, 1975 (1997). 

[6] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, Inc., Orlando, 



11 




FIG. 5: (Color online) The scaled surface tension, 7** 



TcP\ 



2/3 1 



as a function of distance from 



the critical temperature as calculated using the perturbative equation of state displayed as a full 
line (i?* = 00), dashed hue (fi* = 6), dash-dot hue (i?* = 4), dash-dot-dot line (i?* = 3), and 
line-|-circles (i?* = 2.5). 



[7: 

[9 

[10: 
[11 

[12 
[13 
[14 
[15: 

[16 

[17: 



FL, USA, 2001). 

C. H. Bennett, J. Comput. Phys. 22 (1976). 

J. R. Errington and D. A. Kofke, J. Chem. Phys. 127 (2007). 

E. Salomons and M. Mareschal, J. Phys.: Condens. Matter 3 (1991). 

M. J. Haye and C. Bruin, J. Chem. Phys. 100, 556 (1994). 

D. Duque, J. C. Pamies, and L. F. Vega, J. Chem. Phys. 121, 11395 (2004). 
J. J. Potoff and A. Z. Panagiotopoulos, J. Chem. Phys. 112, 6411 (2000). 

J. K. Johnson, J. A. Zollweg, and K. E. Gubbins, Molecular Physics 78, 591 (1993). 

J. A. Barker and D. Henderson, J. Chem. Phys. 47, 4714 (1967). 

J. -P. Hansen and I. McDonald, Theory of Simple Liquids (Academic Press, San Diego, Ca, 

1986). 

E. A. Guggenheim, J. Chem. Phys. 13, 253 (1945). 

J. C. Le Guillou and J. Zinn-Justin, Phys. Rev. Lett. 39, 95 (1977). 



12 



[18] J. C. Le Guillou and J. Zinn-Justin, Phys. Rev. B 21, 3976 (1980). 
[19] H. Okumura and F. Yonezawa, J. Chem. Phys. 113, 9162 (2000). 

[20] J. Perez-Pellitero, P. lingerer, G. Orkoulas, and A. D. Mackie, J. Chem. Phys. 125, 054515 
(2006). 



13 



